C.................... AUR2STR.FOR .............................
C>>>> May 1999: Changes were made to the cascade functions to bring
C>>>> this model into agreement with the FLIP photoelectron model.
C>>>> The altiude dimension was increased to 401 and e* production
C>>>> PRED put in explicitly instead of into cascade production
C>>>> ACASOXT and ACASN2T modified to spread cascade over 3 adjacent bins
C>>>> Excitation potentials in ASETPOT adjusted to bring them in line with
C>>>> values in FLIP 
C.... This is a 2 stream auroral electron model incorporating the
C.... method of Swartz, 1985 (JGR, p6587) for conserving energy
C.... with variable energy grid. See also Stamnes and Rees 1985,
C.... and Banks et al. 1974 (JGR, p1459). This program was first
C.... written by P. Richards in 1982 and modified extensively in
C..... September 1987.
C--N91-  comments like this one refer to modifications in November 1991
C--N91-  The code was modified to include photoelectrons for sunlit auroras
C++++ This program was modified in summer 1997 by Y.Long to allow
C++++ a different way of treating the variable altitude grid. The
C++++ original method is contained in subrotine ATRIS0
C... Version with FPAS conversion put back in ATRIS 9/12/97 
C... I recalled that there is no definition for FPAS and IPAS in 
C... AURORAL PROGRAM, so I just removed that line when I was testing 
C... the program. I've put it back and commented it out, addition to 
C... a little change in keeping the same pattern with ATRISM1.
!... A problem with the inelastic backscatter coefficient was
!... corrected in July 1998. The Solomon values are now used
!... although it is the same for all 3 neutrals. Previously,
!... the coeff for degraded primaries from ionizing collision
!... was was 0.0 and from excitation it was 0.5 (P. Richards)
!... This lowers the LBH peak altitude by about 5 km.
      SUBROUTINE AUR2STR(IZDIM,Z,XN,IMAX,DELKM,ZNE,TE,BG,ZWR,ECHAR
     >    ,ISWFLX,TN,SFAC,EMAX,EMIN,AVMU,DIPANG,ITCONS,TENFLX,SZA)
      PARAMETER (IEDIM=4000,JZDIM=401,IDIM=13)      
      LOGICAL DEBUG
      REAL IONFAC(IEDIM), BM(JZDIM)
      INTEGER IALT(13)
      DIMENSION PEBSC(3),PHIUP(JZDIM),PHIDWN(JZDIM),TSIGNE(JZDIM)
     > ,BG(JZDIM),SIGEX(3),SIGION(3),ZNE(JZDIM),E(IEDIM),TN(IZDIM)
     > ,PRODWN(IEDIM,JZDIM),Z(JZDIM),PRODUP(IEDIM,JZDIM),SIGEL(3)
     > ,XN(3,JZDIM),T1(JZDIM),TE(JZDIM),T2(JZDIM),A(JZDIM)
     > ,B(JZDIM),C(JZDIM),D(JZDIM),SECPRD(3,JZDIM),ENABEX(3,JZDIM)
     > ,ENABIO(3,JZDIM),PROBSE(IEDIM),EHEAT(JZDIM),PHISUM(JZDIM)
     > ,DELTE(IEDIM),IE(IEDIM),JDP(3,IEDIM),JIE(IEDIM)
     > ,IEDUM(IEDIM),EXIFAC(3,IEDIM),ALPDUM(JZDIM),GAMDUM(JZDIM)
     > ,SUMEX(3),SUMIO(3),EPOTIN(3,IEDIM),ELOSSX(3,IEDIM)
     > ,ELPRIM(IEDIM),SPRD(3,6),PHIN(IEDIM),ENERAL(JZDIM),HE(JZDIM)
     > ,SUMION(JZDIM)
      REAL DS(JZDIM), PRED(JZDIM),EM_LBH(JZDIM),EM_O1D(JZDIM)
      REAL ETIROS(IDIM),PHITIROS(IDIM),RIONTOT(JZDIM)   !.. Arrays for TIROS factors
      
        DATA  JMAX,M ,M2,IWR,M1,SHAPE,RECRPI, DEBUG,     PI
     >    /     0 ,41,2 , 1 ,1 , 14.0, 0.564, .FALSE.,3.14159  /
      !--- Branching ratios for ion states were updated Sep 91 by P. Richards
      !--- N2+ from Doering and Goembel, JGR 1991, page 16025. fraction of
      !--- N+ from Richards and Torr JGR 1985, page 9917. fraction of O+ from
      !--- O2 dissociation is from Rapp and Englander Golden J. Chem Phys, 1965
      !--- We still need to find cross sections for the higher levels of O+
      DATA SPRD/.4,.47,.42, .4,.28,.34, .2,.15,.08,0.0,.05,0.0,0.0,
     >   .05,0.0,0.0,0.0,0.16/
      DATA BM/JZDIM*1.0/      ! not used in this program.
      DATA PRED/JZDIM*0.0/    ! not used in this program.
C
      !.. check to see if dimensions are OK
      IF(JZDIM.NE.IZDIM) THEN
         WRITE(6,185)
         CALL RUN_ERROR    !.. print ERROR warning in output files
         STOP
      ENDIF
 185  FORMAT(/2X,' ** ERROR:JZMAX must equal IZMAX in routine AUR2STR'/)

      !.. dip angle factor, and altitude step in cm
      DIPMU=AVMU*SIN(DIPANG)
      DELZ=DELKM*1.0E+5
c      DO I = M2, IMAX-1
      DO I = M2, IMAX
        DS(I)=DELZ/BG(I)
        !..DS(I)=DELZ/BG(I)
      ENDDO
C
       CALL ALFLUX(IMAX,JZDIM,Z,INALT,IALT)

      !...... Setting up energy cell boundaries first time thru only
         CALL ACELLS(IEDIM,EMAX,SHAPE,JMAX,E,DELTE)
C
      !... set up bins for degraded primaries from ionizations, excitations
         DO 3 J=1,JMAX
      !.. Find the ionization and excitation energies to find energy loss
      !.. Note that only EPOTIN(3,J) for N2 is actually used
         CALL ASETPOT(IEDIM,J,E(J),EPOTIN,ELOSSX)
      !.. ASETPRI calculates the probabilities for the secondary electrons
      !.. from ionizing collisions. Here it is just used to return the
      !.. average energy of the secondaries (AVESEC) to calculate
      !.. the energy of the degraded primaries (ELPRIM).
         CALL ASETPRI(IEDIM,J,JMAX,E,DELTE,EPOTIN(3,J),IEDUM,PROBSE,ISEC
     > ,AVESEC)
C
         ELPRIM(J)=EPOTIN(3,J)+AVESEC
      !... allocate bins (IE(J)) for degraded primaries for ionization
      !... and calculate Swartz factor IONFAC(J). If the degraded 
      !... electrons don't go in the next lowest bin, they are spread
      !... over 2 bins to smooth out the spectrum. See also ACASION
         IONFAC(J)=1.0
         CALL AFNDBIN(J,JMAX,ELPRIM(J),E,IE)
         if (ie(j).ne.0) then      !xiang added

            EK=E(IE(J))
            IF(J.LT.JMAX) EK=(2.0*E(IE(J))+E(IE(J)+1))/3.0
            IF(IE(J)-J.GT.1) EK=(2*E(IE(J))+E(IE(J)+1)+3*E(IE(J)-1))/6.0
            IF(IE(J).GT.0.AND.IE(J).LE.JMAX)
     >         IONFAC(J)=ELPRIM(J)/(E(J)-EK)
         endif                     !xiang added
      !... allocate bins for degraded primaries for excitation of O and N2
      !... JIE(J), and calculate Swartz factor EXIFAC(J)
      DO 5 I=1,3
         EXIFAC(I,J)=1.0
         CALL AFNDBIN(J,JMAX,ELOSSX(I,J),E,JIE)
         JDP(I,J)=JIE(J)
         IF(JIE(J).GT.0.AND.JIE(J).LE.JMAX)
     >      EXIFAC(I,J)=ELOSSX(I,J)/(E(J)-E(JDP(I,J)))
 5    CONTINUE
      !...... diagnostic printout
         IF(DEBUG.AND.(J.EQ.JMAX.OR.MOD(J,15).EQ.0)) WRITE(104,200)
 200     FORMAT(2X,'J   IE JDP JDP JDP     E      DELTE    IONFAC ',
     >  '  EXIFAC   EXIFAC     E        E        E')
         IF(DEBUG) WRITE(104,195) J,IE(J),(JDP(I,J),I=1,3),E(J),DELTE(J)
     >    ,IONFAC(J)
     >    ,EXIFAC(1,J),EXIFAC(2,J),E(IE(J)),E(JDP(1,J)),E(JDP(3,J))
         IF(DEBUG.AND.(J.EQ.JMAX.OR.MOD(J,15).EQ.0)) WRITE(154,207)
 207     FORMAT(2X,'J   IE JDP JDP JDP     E      DELTE',
     >  ' ELPRIM , EPOTIN , EPOTIN , ELOSSX , ELOSSX , AVESEC')
         IF(DEBUG) WRITE(154,195) J,IE(J),(JDP(I,J),I=1,3),E(J),DELTE(J)
     > ,ELPRIM(J),EPOTIN(1,J),EPOTIN(3,J),ELOSSX(1,J),ELOSSX(3,J)
     >    ,AVESEC
 195     FORMAT(5I4,1P,22E9.1)
 3    CONTINUE

      !.. CALL ADIAGS(IEDIM,JMAX,ECHAR,E,EPOTIN,ELPRIM,ELOSSX)
      
      PHINSUM=0.0
      ESUM=0.0
      EBSC=0.0
      PHIBSC=0.0
      IWR=1
C
      !..... initialization of arrays and find index for print altitude IWR
      DO 13 I=M1,IMAX
        IF(Z(I).LE.ZWR) IWR=I
        EHEAT(I)=0.0
        EM_LBH(I)=0.0
        EM_O1D(I)=0.0
        PHIUP(I)=0.0
        PHIDWN(I)=0.0
        SUMION(I)=0
        DO 12 IS=1,3
          ENABEX(IS,I)=0.0
          ENABIO(IS,I)=0.0
          SECPRD(IS,I)=0.0
 12     CONTINUE
 13   CONTINUE

      !... Used to test energy conservation in local equilibrium
      !... ITCONS=0 - normal aurora; ITCONS=1 - local production
      !... ITCONS=2 - local production + local equilibrium
      DO 21 JJ=1,JMAX
      DO 21 I=M1,IMAX
         IF(ITCONS.NE.0) THEN
             PRODUP(JJ,I)=1.0
             PRODWN(JJ,I)=1.0
         ELSE
             PRODUP(JJ,I)=0.0
             PRODWN(JJ,I)=0.0
         ENDIF
 21   CONTINUE

      !... write altitudes for auroral fluxes
      IF(ZWR.GT.80.AND.ZWR.GE.Z(1)) THEN
         WRITE(6,317) Z(IWR)
      ENDIF
 317   FORMAT(/5X,'Auroral flux at',F8.2,'km altitude')

      !................ Main energy loop begins ........................
C
      DO 82 J=1,JMAX
      !IF(E(J).LT.1000) GO TO 80
C
      !.... Get elastic cross sections from ELASTC and inelastic from SIGEXS
      CALL ELASTC(E(J),PEBSC,SIGEL)
      CALL SIGEXS(E(J),SIGEX,SIGION,SIGO1D)
      CALL A1PIXS(E(J),SIGLBH)               !... LBH cross section for test
      CALL BACKSCAT(E(J),BSCAT)              !... inelastic backscatter coeff

      !.... set up the probability distribution for the secondary electrons
      !.... (PROBSE) and find energy bins for the highest energy secondary
      !.... (ISEC) and degraded primaries 
      CALL ASETPRI(IEDIM,J,JMAX,E,DELTE,EPOTIN(3,J),IEDUM,PROBSE,ISEC
     > ,ESEC)

      !.... Energy loss to coulomb collisions with thermal electrons
      DO 11 I=M1,IMAX
      ET=8.618E-5*TE(I)
      TSIGNE(I)=((3.37E-12*ZNE(I)**0.97)/(E(J)**0.94))*((E(J)-ET)/
     >   (E(J)-(0.53*ET)))**2.36/(E(J)-E(J+1))
      TEST=TSIGNE(I)*DS(I)
      !.... check FLIP      TSIGNE(I)=TSIGNE(I)/DELTE(J)
      IF(TEST.GT.1.0)TSIGNE(I)=1.0/(E(J)-E(J+1))/(DS(I))
 11   CONTINUE
C
C.... For T1 and T2 see Nagy and Banks 1970 or B&K P258 ............
      DO 30 I=M1,IMAX
      T1(I)=0.0
      T2(I)=TSIGNE(I)
      DO 31 IS=1,3
      T1(I)=T1(I)+XN(IS,I)*SIGEL(IS)*PEBSC(IS) 
      TABSXS=SIGEX(IS)*EXIFAC(IS,J)+SIGION(IS)*IONFAC(J)
 31   T2(I)=T2(I)+XN(IS,I)*(SIGEL(IS)*PEBSC(IS)+TABSXS)    
C...      WRITE(103,311) I,Z(I),T1(I),T2(I),TSIGNE(I)      
 30    CONTINUE
C
C
      !..... Set boundary conditions on fluxes ......
      M=IMAX-1
      !..... lower boundary
       PHIDWN(M1)=0.0
       PHIUP(M1)=0.0
C
C ----------------  Incident energy flux --------------------
      !..The Gaussian spectrum equation from Daniel and Strickland 1986.
      !.. TENFLX = QM or QG, ECHAR = Em or Eg, WEQ = W
      IF(SFAC.GT.0.01) THEN
         WEQ=SFAC*ECHAR
         PHIN(J)=0.564*TENFLX*
     >      EXP(-(E(J)-ECHAR)**2/WEQ**2) / (WEQ*ECHAR)
      !..The Maxwellian spectrum equation from Daniel and Strickland 1986.
      !.. Note that the 0.5 factor comes from the integral of the distribution 
      !.. over energy = E**2*EXP(-E/E0)dE= 
      !..          -E**2*E0*EXP(-E/E0)+2.0*E0**3*EXP(-E/E0)[-E/E0-1]
      !.. Taking the limit from 0 to infinity the integral = 2*E0**3
      !.. See Beyer Standard Math Tables for the integral of E**2*EXP(-E/E0)
      ELSE 
         PHIN(J)= 0.5 *TENFLX* E(J)*EXP(-E(J)/ECHAR)/ECHAR**3
      ENDIF

      !.. CALL EVANS_FLUX(E(J),PHIN(J))        !.. Evans Flux
      !.. CALL FAST_FLUX(E(J),PHIN(J))         !.. FAST Flux
      !.. CALL TJFRE_FLUX(E(J),PHIN(J))        !.. Fuller-Rowell and Evans Flux
      !.. CALL TIROHIST(IR,E(J),AVMU,PHIN(J))  !.. TIROS flux from F-R&E 1987 Fig.1
      !.. PHINREES puts fluxes into TIROS bins for REES routine. Modify TIROHIS
      !.. to select a single bin
       CALL PHINREES(E(J),DELTE(J),PHIN(J),PHITIROS,ETIROS)

       !.. Low energy tail flux from Meier et al 1989
       BFAC=1/(0.8*ECHAR)
       IF(ECHAR.GE.500) BFAC=1/(0.1*ECHAR+0.35)
       TAILFLUX=0.4 * 0.37* TENFLX*EXP(-BFAC*E(J))/E(J)/ECHAR
       TAILFLUX=0.1*TENFLX/ECHAR/E(J)
       TAILFLUX=0.001*TENFLX*(1/E(J)/E(J)+4E-7*SQRT(E(J)))

       ESUM=ESUM+E(J)*PHIN(J)*DELTE(J)      !.. Total Energy precipitating
       PHINSUM=PHINSUM+PHIN(J)*DELTE(J)     !.. Total Flux precipitating

C --------------------- END incident fluX ----------------------
       IF(ITCONS.EQ.0) PHIDWN(IMAX)=PHIN(J)/AVMU
     >    +1.0*PHIUP(M-1) !.. add BS flux
c     > +TAILFLUX

C
C.... Pitch is taken from FLIP it changes the pitch angle, but also
C.... normalizes T1, T2 etc. The pitch angle does not vary at low
C.... altitudes
       CALL APITCH(J,0,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z,DIPANG
     >           ,AVMU)
C
      !.. setting up the coefficients of the parabolic PDE for down flux.
       CALL ATRIS1(1,M2,M,J,M1,BM,BG,Z,IMAX+1,DELZ,PRED,PRODUP,PRODWN,
     &  PHIDWN,PHIUP,T1,T2,DS)

       CALL APITCH(J,1,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z,DIPANG
     >           ,AVMU)
      !... local equilibrium fluxes
       IF(ITCONS.EQ.2) THEN
          DO 17 I=1,IMAX
          PHIUP(I)=PRODUP(J,I)/(T2(I)-T1(I))
          PHIDWN(I)=PRODWN(J,I)/(T2(I)-T1(I))
 17       CONTINUE
      !..... sum the energy in the backscattered flux
       ELSE
          EBSC=EBSC+E(J)*PHIUP(IMAX)*DELTE(J)*AVMU
          PHIBSC=PHIBSC+PHIUP(IMAX)*DELTE(J)*AVMU
       ENDIF
C
      !... This loop for calculating energy absorption etc.
      IFYNEG=0
      DO I=M1,IMAX
        IF(PHIUP(I).LT.0.0.OR.PHIDWN(I).LT.0.0) IFYNEG=I
        PHISUM(I)=(PHIUP(I)+PHIDWN(I))*DELTE(J)
        !.. sum energy going into ions (ENABIO) and excitations (ENABEX)
        DO 130 IS=1,3
        IF(E(J).GT.ELOSSX(IS,J)) ENABEX(IS,I)=ENABEX(IS,I)+
     >    ELOSSX(IS,J)*PHISUM(I)*SIGEX(IS)*XN(IS,I)
        IF(E(J).GT.EPOTIN(IS,J)) ENABIO(IS,I)=ENABIO(IS,I)+
     >    EPOTIN(IS,J)*PHISUM(I)*SIGION(IS)*XN(IS,I)
        !.. running total of energy absorbed
        IF(I.EQ.IWR) ERUN=ERUN+EPOTIN(IS,J)*PHISUM(I)*SIGION(IS)*
     >    XN(IS,I)*EXIFAC(IS,I)+ELOSSX(IS,J)*PHISUM(I)*SIGEX(IS)*
     >    XN(IS,I)*IONFAC(IS)
        !... SECPRD is secondary ion production rate
        SUMION(I)=SUMION(I)+PHISUM(I)*SIGION(IS)*XN(IS,I)
 130    SECPRD(IS,I)=SECPRD(IS,I)+PHISUM(I)*SIGION(IS)*XN(IS,I)
        !..... Thermal electron heating rate
        EHEAT(I)=EHEAT(I)+PHISUM(I)*TSIGNE(I)
        IF(DEBUG.AND.I.EQ.IWR) THEN
          ERUN=ERUN+PHISUM(I)*TSIGNE(I)
          IF(DEBUG) WRITE(104,311) I,Z(I),ERUN,EHEAT(I)
     >     ,(ENABEX(IS,I),IS=1,3),(ENABIO(IS,I),IS=1,3)
        ENDIF
      ENDDO

        IF(IFYNEG.NE.0) WRITE(6,196) NINT(Z(IFYNEG))
 196  FORMAT('  ** Negative flux at altitude= ',I6,' km')

      !.... Loop to calculate cascade production
      DO 83 I=M1,IMAX

      !.... calculate  cascade production from  degraded primaries in
      !.... ionizing collisions
      IF(E(J).GT.EPOTIN(3,J)) THEN
      CALL ACASION(JZDIM,IEDIM,J,I,IWR,Z(I),XN,PRODUP,PHIUP(I),PHIDWN(I)
     > ,E,DELTE,SIGION,JMAX,IE(J),IONFAC(J),BSCAT)
      CALL ACASION(JZDIM,IEDIM,J,I,IWR,Z(I),XN,PRODWN,PHIDWN(I),PHIUP(I)
     > ,E,DELTE,SIGION,JMAX,IE(J),IONFAC(J),BSCAT)
      !.. distribute the secondaries over low energy bins
      CALL ACASSEC(JZDIM,IEDIM,J,I,IWR,Z(I),XN,PRODUP,PHIUP(I),PHIDWN(I)
     > ,E,DELTE,SIGION,PROBSE,ISEC,JMAX,0.5)
      CALL ACASSEC(JZDIM,IEDIM,J,I,IWR,Z(I),XN,PRODWN,PHIDWN(I),PHIUP(I)
     > ,E,DELTE,SIGION,PROBSE,ISEC,JMAX,0.5)
      ENDIF
C
      !..... Cascade from thermal electron collisions. DELTE/DELTE
      !..... takes into account differing energy cells
      TSIGNE(I)=TSIGNE(I)*DELTE(J)/DELTE(J+1)
      DATA EBSCX/0.0/
      PRODUP(J+1,I)=PRODUP(J+1,I)+PHIUP(I)*TSIGNE(I)*(1-EBSCX)
     > +PHIDWN(I)*TSIGNE(I)*EBSCX
      PRODWN(J+1,I)=PRODWN(J+1,I)+PHIDWN(I)*TSIGNE(I)*(1-EBSCX)
     > +PHIUP(I)*TSIGNE(I)*EBSCX
C
      !..... calculate cascade from atomic O
      CALL ACASOXT(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODUP,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JDP(1,J),EXIFAC(1,J),BSCAT)
      CALL ACASOXT(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODWN,PHIDWN
     > ,PHIUP,E,DELTE,SIGO1D,SIGION,JMAX,JDP(1,J),EXIFAC(1,J),BSCAT)
      !... cascade from N2
      CALL ACASN2T(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODUP,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JDP(3,J),EXIFAC(3,J),BSCAT)
      CALL ACASN2T(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODWN,PHIDWN
     > ,PHIUP,E,DELTE,SIGO1D,SIGION,JMAX,JDP(3,J),EXIFAC(3,J),BSCAT)
      !... cascade from O2
      CALL ACASO2T(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODUP,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JDP(2,J),EXIFAC(2,J),BSCAT)
      CALL ACASO2T(JZDIM,IEDIM,J,I,IWR,Z(I),SIGEX,XN,PRODWN,PHIDWN
     > ,PHIUP,E,DELTE,SIGO1D,SIGION,JMAX,JDP(2,J),EXIFAC(2,J),BSCAT)


      !... emission rates for diagnostics
      EM_LBH(I)=EM_LBH(I)+(PHIUP(I)+PHIDWN(I))*SIGLBH*XN(3,I)
      EM_O1D(I)=EM_O1D(I)+(PHIUP(I)+PHIDWN(I))*SIGO1D*XN(1,I)
 83   CONTINUE

      !.. IF(ZWR.GT.80) WRITE(6,313) E(J),(PHIUP(K),K=IWR-3,IWR+3),
      !..>        (PHIDWN(K),K=IWR-3,IWR+3)

      !......... Diagnostic Outputs............
      IF(Z(IWR).GT.80) THEN
        IF(J.EQ.1) WRITE(6,181) NINT(Z(IWR))
 181    FORMAT(/2X,' J, Energy, input flux, up flux,'
     >   ,' down flux, up and down production rates at altitude='
     >   ,I7,' km')
        IF(J.EQ.1.OR.MOD(J,15).EQ.0) WRITE(6,193)
 193    FORMAT(5X,'J      EV      PRECIP     PHIUP    PHIDWN    PRODUP'
     > ,'   PRODDWN    SIGOX      SIGO2     SIGN2    ISIGOX    ISIGO2'
     > ,'    ISIGN2    PEBSOX    PEBSO2    PEBSN2    BSCAT     SIGEOX'
     > ,'    SIGEO2    SIGEN2   TAILFLUX   TPREC    SUMION     EHEAT')
        WRITE(6,95) J,E(J),PHIN(J)/AVMU,PHIUP(IWR),PHIDWN(IWR)
     >   ,PRODUP(J,IWR),PRODWN(J,IWR),(SIGEX(IS),IS=1,3)
     >   ,(SIGION(IS),IS=1,3),(PEBSC(IS),IS=1,3),BSCAT,
     >   (SIGEL(IS),IS=1,3),TAILFLUX,(PHIN(J)/AVMU+TAILFLUX),SUMION(IWR)
     >   ,EHEAT(IWR)
      ENDIF
 
      IF(DEBUG) THEN
        WRITE(135,395) E(J),(PHIUP(IALT(L)),L=1,INALT)
        WRITE(136,395) E(J),(PHIDWN(IALT(L)),L=1,INALT)
395     FORMAT(X,F7.1,1P22E9.2)
        !..........write  cross sections
        IF(J.EQ.1) WRITE(113,184)
 184    FORMAT(//2X,' Excitation and ionization cross sections '/)
        IF(J.EQ.1.OR.MOD(J,15).EQ.0) WRITE(113,192)
 192    FORMAT(4X,'J    E(J)   SIGOX    SIGO2    SIGN2    ISIGOX',
     >   '   ISIGO2   ISIGN2  SIGELOX  SIGELO2  SIGELN2   SIGO1D'
     >   ,1X,' EL_BSC   INEBSC')
        WRITE(113,311) J,E(J),(SIGEX(IS),IS=1,3),(SIGION(IS),IS=1,3)
     >  ,(SIGEL(IS),IS=1,3),SIGO1D,PEBSC(1),BSCAT
      ENDIF

      !.... PRODWN is used to store FLUX*DELTE for output at the end
      DO 902 I=M1,IMAX
 902  PRODWN(J,I)=PHISUM(I)
 82   CONTINUE
 80   CONTINUE
      !.................. END OF ENERGY LOOP
C
C
      !.... diagnostic printout of altitude profiles
      IF(DEBUG) THEN
         DO 87, I=M1+1,IMAX-1
           WRITE(144,'(F8.2,1P,9E10.2)') Z(I),EM_O1D(I),EM_LBH(I)
     >        ,XN(1,I),XN(3,I),SECPRD(1,I),SECPRD(2,I),SECPRD(3,I)
 87      CONTINUE
      ENDIF
      !...... IF EMIN > 0 write out low energy secondary production
      IF(EMIN.GT.0) THEN
         WRITE(6,*) '  J= ',J
         WRITE(6,197) (PRODWN(IJ,IWR),IJ=JMAX,JMAX-10,-1)
 197     FORMAT(/2X,' PRODWN(JMAX)-> PRODWN(JMAX-10)',/1P,22E9.1/)
         WRITE(6,198) (PRODWN(IJ,IWR),IJ=1,10)
 198     FORMAT(/2X,' PRODWN(1)-> PRODWN(10)',/1P,22E9.1/)
      ENDIF

      !.. This calculate and outputs the energy balance for auroral precipitation
      IF(ZWR.GT.80) THEN
        CALL CONSERVE(M1,IMAX,JZDIM,ZWR,Z,ECHAR,TENFLX,EBSC,ESUM,EHEAT,
     >    ENABEX,ENABIO,SECPRD)
        WRITE(6,*) '     ZALT     RHO       ZMASS     Ntot       Z'
     >   ,'       R         ZoR      LAMBDA   REESION    TOTION'
     >   ,'    SUMION   PHINSUM    IONRAT    RIONTOT'

        IR=0
        IF(IR.GT.0) THEN     !.. Section begin for REES calculation          
          DO IR=1,IDIM
            RREES=4.57E-6*(ETIROS(IR)/1000)**1.75
            DO I=M1+1,IMAX-1
              ZMASS=1.67E-24*(16*XN(1,I)+32*XN(2,I)+28*XN(3,I))/
     >          (XN(1,I)+XN(2,I)+XN(3,I))
              GRAV=980.665/(1.0+Z(I)/6370.0)**2  !.. gravity
              SHGT=1.38E-16*TN(I)/GRAV/ZMASS
              RMASS=(XN(1,I)+XN(2,I)+XN(3,I))*ZMASS*SHGT
              IF(RMASS.GT.RREES) THEN
                IM=I
                RM=RMASS
              ENDIF
            ENDDO
            RRHO=1.67E-24*(16*XN(1,IM)+32*XN(2,IM)+28*XN(3,IM))
            RDEN=(XN(1,IM)+XN(2,IM)+XN(3,IM))
      
            DO I=M1+1,IMAX-1
             CALL REESAUR(Z(I),TN(I),XN(1,I),XN(2,I),XN(3,I),ETIROS(IR),
     >       TENFLX,PHITIROS(IR),RRHO,RDEN,SUMION(I),REESION,RIONTOT(I))
             RIONTOT(I)=RIONTOT(I)+REESION*PHITIROS(IR)
            ENDDO
          ENDDO
        ENDIF     !.. Section end for REES calculation 
      
      ENDIF    ! End ZWR
      
      ENSUM=0.0
      !----- calculate the total energy going into each altitude for FLIP
      DO I=M1+1,IMAX-1
         ENERAL(I)=EHEAT(I)+ENABEX(1,I)+ENABEX(2,I)+ENABEX(3,I )
     >     +ENABIO(1,I)+ENABIO(2,I)+ENABIO(3,I)
         !..WRITE(6,95) I,Z(I),ENERAL(I),EHEAT(I),ENABEX(3,I),ENABIO(3,I)
         ENSUM=ENSUM+SQRT(ENERAL(I)*ENERAL(I))*1.0E5*(Z(I+1)-Z(I))
      !..  WRITE(93,'(F10.2,1P,9E10.2)')  Z(I),ENERAL(I),ENSUM/6.25E11
      ENDDO
      IF(ZWR.GT.80) WRITE(6,'(A,1PE10.2)') ' Total E=',ENSUM/6.25E11
    
      !.. write auroral production frequencies for FLIP
      CALL AURFRE(IMAX,JMAX,JZDIM,IEDIM,Z,E,PRODWN,ENERAL)

 95   FORMAT(1X,I6,F10.1,1P29E10.2)
 171  FORMAT(X,F8.1,1P22E11.2)
 311  FORMAT(1H ,I5,F7.0,1P22E9.1)
 313  FORMAT(F9.1,1P22E9.2)
          END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::      
        SUBROUTINE ATRIDAG(JZDIM,DELTA,IF,L,A,B,C,D,ALPHA,GAMMA)
C.... FOR SOLVING A SYSTEM OF LINEAR SIMULTANEOUS EQUATIONS WITH A       
C     TRIDIAGONAL COEFF MATRIX. THE EQNS ARE NUMBERED FROM IF TO L, 
C     & THEIR  SUB-DIAG. , DAG. ,& SUPER-DIAG COEFFS. ARE STORED
      !..    IN THE ARRAYS A ,B& C
       DIMENSION A(JZDIM),B(JZDIM),C(JZDIM),D(JZDIM)
       DIMENSION ALPHA(JZDIM),DELTA(JZDIM),GAMMA(JZDIM)
      !..  COMPUTE INTERMEDIATE ARRAYS ALPHA & GAMMA
        ALPHA(IF)=B(IF)
        GAMMA(IF)=D(IF)/ALPHA(IF)
        IFP1=IF+1
        DO 1 I=IFP1,L
        ALPHA(I)=B(I)-A(I)*C(I-1)/ALPHA(I-1)
        GAMMA(I)=(D(I)-A(I)*GAMMA(I-1))/ ALPHA(I)
1       CONTINUE
      !..  COMPUTE FINAL SOLUTION VECTOR V
        DELTA(L)=GAMMA(L)
        LAST=L-IF
        DO 2 K=1,LAST
        I=L-K
        DELTA(I)=GAMMA(I)-C(I)*DELTA(I+1)/ALPHA(I)
2       CONTINUE
        RETURN
        END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C......... Subroutine to set energy grid for the 2 stream program
C......... EMAX is the maximum energy required. JMAX= number of energy cells
C......... It returns the midpoint energies E and the width of the energy
C......... cells DELTE. 
      SUBROUTINE ACELLS(IEDIM,EMAX,SHAPE,JMAX,E,DELTE)
      DIMENSION E(IEDIM),DELTE(IEDIM)
C
       !... If EMAX is negative use 1eV steps but check dimension first
       IF(EMAX.LE.-IEDIM) THEN
          WRITE(6,90) ABS(EMAX),IEDIM
          CALL RUN_ERROR    !.. print ERROR warning in output files
          STOP
       ENDIF
 90      FORMAT(/1X,' *** Cannot solve because ABS(EVMAX)='1P,E10.1
     >          ,' is bigger than energy dimension=',I6)

      !...... set up cell midpoint energies according to energy lost
      !.... in ionization
      E(1)=ABS(EMAX)+0.5
         DO 75 J=1,IEDIM
C
      !.. If EMAX < 101 use 1 eV energy steps else large steps
      IF(EMAX.LT.101) THEN
         E(J+1)=E(J)-1.0
         DELTE(J)=1.0
      ELSE
C
      !.....  LINEAR energy grid between energies E1 and E2
      !..... where the energy steps are S1 and S2
         DATA E1,E2,S1,S2/20.0,20000,1.0,500.0/
      !... find gradient and abscissa, then DELTE
         GRAD=(S2-S1)/(E2-E1)
         ABSC=S1-E1*GRAD
         DELTE(J)=(GRAD*E(J)+ABSC)
         IF(DELTE(J).LT.1.1) THEN
            GRAD=0.0
            ABSC=1.0
            DELTE(J)=1.0
         ENDIF
         E(J+1)=(E(J)-0.5*DELTE(J)-0.5*ABSC)/(1+0.5*GRAD)
         IF(E(J+1).LT.0.0) GO TO 76
         E(J+1)=INT(E(J+1))+0.5
      ENDIF
         !..WRITE(6,95) J,E(J),E(J+1),DELTE(J)
 95      FORMAT(I5,9F9.1)
 75      CONTINUE
 76      CONTINUE
         JMAX=J
C
      EMAX=ABS(EMAX)       !... make sure EMAX is positive
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C........ This is the program for calculating the probability (PROB) per eV
C........ for a secondary electron to have an energy ELECEN
      SUBROUTINE ASETPRI(IEDIM,J,JMAX,E,DELTE,EPOTI,IE,PROB,ISEC,ESEC)
      DIMENSION E(IEDIM),DELTE(IEDIM),IE(IEDIM),PROB(IEDIM)
      ESEC=0.0
      ISEC=0
      DO 75 IK=1,IEDIM/2
 75   PROB(IK)=0.0
      IF(E(J).LE.EPOTI+1.0) RETURN
      !.... energy lost by lowest energy primary
      ENLSS=0.5*(E(J)+EPOTI)
      !.. Find max secondary energy bin. 
      CALL AFNDBIN(J,JMAX,ENLSS,E,IE)
      ISEC=IE(J)
      IF(ISEC.LT.1.OR.ISEC.GT.JMAX) RETURN
C
      !.. Probability of secondary electron of energy E(IK)
      !.. See Banks et al. 1974
        PROBS=0.0
        ESEC=0.0
        DO 115 IK=ISEC,JMAX
        PROB(IK)=1.0/(1+(E(IK)/14.0)**2)/ATAN(E(ISEC)/14.0)/14.0
        ESEC=ESEC+E(IK)*PROB(IK)*DELTE(IK)
        PROBS=PROBS+PROB(IK)*DELTE(IK)
      !..        WRITE(6,90) IK,E(IK),DELTE(IK),PROB(IK),PROBS,ESEC
C.. 90   FORMAT(I6,1P,9E10.2)
 115    CONTINUE
      !.. total energy in secondary electrons (ESEC)
        ESEC=ESEC/PROBS
      !.. divide by total probability (PROBS) to ensure it =1
      IF(PROBS.LE.0.0) RETURN
      DO 215 IK=ISEC,JMAX
 215  PROB(IK)=PROB(IK)/PROBS
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE AFNDBIN(J,JMAX,ELOSS,EE,IE)
C........ A program to determine which bin to allocate a degraded electron
C........ J= index of the primary of energy EE(J). JMAX is the max J, ELOSS=
C........ energy lost by primary and the index is returned in IE
      INTEGER J,JMAX,IE(JMAX)
      REAL EE(JMAX)
      !.. Energy of degraded primary
      EDPRIM=EE(J)-ELOSS
      !..... Test to see if energy loss is too great
      IF(EDPRIM.LT.EE(JMAX)) THEN
        IE(J)=0
      ELSE
      !.... Test degraded energy against next lowest bin until it fits
      DO 29 I=J,JMAX-1
      IF(EDPRIM.GT.EE(I)) GO TO 30
      IF(ABS(EDPRIM-EE(I)).LE.ABS(0.5*(EE(I)-EE(I+1)))) GO TO 30
C...      WRITE(6,90) I,EDPRIM,EE(I)
C... 90   FORMAT(I5,9F9.1)
 29   CONTINUE
 30   CONTINUE
      !.... store index for return
      IE(J)=I
      IF(IE(J).EQ.J) IE(J)=J+1
      !..      WRITE(6,90) IE(J),EE(IE(J))
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C......... This program returns the average ionization potential and
C......... average excitation energy for O, O2, and N2 for an electron
C.......... of energy ELECEN. The energy loss for O excitation does not
C.......... include O(1D) which is treated separately. Similarly, N2
C.......... does not include N2*
      SUBROUTINE ASETPOT(IEDIM,J,ELECEN,EPOTIN,ELOSSX)
      DIMENSION EPOTIN(3,IEDIM),ELOSSX(3,IEDIM)
C
      !........ Average excitation potentials. First O
      IF(ELECEN.LE.10) THEN
         ELOSSX(1,J)=2.0
      ELSE
         ELOSSX(1,J)=2.0+10.0*SQRT(1.0-10.0/ELECEN)   ! May 1999
      ENDIF
      !.. N2
      IF(ELECEN.LE.6.0) THEN
         ELOSSX(3,J)=1.0
      ELSE 
          ELOSSX(3,J)=4.0+7.0*SQRT(1.0-6.0/ELECEN)   ! May 1999
      ENDIF
      !..... set O2=N2
      ELOSSX(2,J)=ELOSSX(3,J)
      !... ionization potentials for O, O2, N2 are equal to that
      !... for N2. O should maximize at about 18 eV instead of 20eV
      IF(ELECEN.GT.15.0) EPOTIN(1,J)=10+10.0*SQRT(1-15.0/ELECEN)
      IF(ELECEN.GE.13.6.AND.EPOTIN(1,J).LT.13.6) EPOTIN(1,J)=13.6
      IF(EPOTIN(1,J).GT.18) EPOTIN(1,J)=18.0
      IF(ELECEN.GE.15) EPOTIN(2,J)=10+10.0*SQRT(1-15.0/ELECEN)
      IF(ELECEN.GE.12.6.AND.EPOTIN(2,J).LT.12.6) EPOTIN(2,J)=12.6
      IF(ELECEN.GE.15.0) EPOTIN(3,J)=10+10.0*SQRT(1-15.0/ELECEN)
      IF(ELECEN.GE.15.6.AND.EPOTIN(3,J).LT.15.6) EPOTIN(3,J)=15.6
      IF(ELECEN.LT.13.6) EPOTIN(1,J)=0
      IF(ELECEN.LT.12.6) EPOTIN(2,J)=0
      IF(ELECEN.LT.15.6) EPOTIN(3,J)=0
      RETURN
      END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C...... to calculate cascade production rates from O - O(1d)
C...... JOE= index of bin for depositing degraded electrons
      SUBROUTINE ACASOXT(JZDIM,IEDIM,J,IZ,IWR,ALT,SIGEX,XN,PROCAS,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JOE,EXIFAC,BSCX)
      DIMENSION SIGEX(3),XN(3,JZDIM),PROCAS(IEDIM,JZDIM),PHIUP(JZDIM)
     > ,E(IEDIM),DELTE(IEDIM),PHIDWN(JZDIM),SIGION(3)
C
      !...... cascade from O - O(1D)
      IF(JOE.GT.0.AND.JOE.LE.JMAX.AND.E(J).GE.10.0)  THEN
         JOEM1=JOE
         JOEM2=JOE
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-1) JOEM1=JOE+1  !.. May 1999
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-2) JOEM2=JOE+2  !.. May 1999
         CORSIG=(SIGEX(1)-SIGO1D)*EXIFAC
         PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)
      !.. allocate production from bin J in 2 adjacent bins if small bins
         PREX=0.33333*PREX*CORSIG*XN(1,IZ)
         PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+PREX/DELTE(JOE)
         PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+PREX/DELTE(JOEM1)
         PROCAS(JOEM2,IZ)=PROCAS(JOEM2,IZ)+PREX/DELTE(JOEM2)
      ENDIF
C
      !...... cascade from O(1D) 
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)
      JK=J+2
      IF(JK.LE.JMAX) THEN
         IF(DELTE(J).GT.1) JK=J+1
         CORSIG=2.0/(E(J)-E(JK))
         PROCAS(JK,IZ)=PROCAS(JK,IZ)+PREX*SIGO1D*CORSIG*XN(1,IZ)/
     >   DELTE(JK)
      ENDIF
      IF(J.EQ.-979) WRITE(6,*) IWR  !.. Dummy statement to use IWR
        RETURN
        END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C........ to calculate cascade production rates from N2
C........ JOE= index of bin for depositing degraded electrons
      SUBROUTINE ACASN2T(JZDIM,IEDIM,J,IZ,IWR,ALT,SIGEX,XN,PROCAS,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JOE,EXIFAC,BSCX)
      DIMENSION SIGEX(3),XN(3,JZDIM),PROCAS(IEDIM,JZDIM),PHIUP(JZDIM)
     > ,E(IEDIM),DELTE(IEDIM),PHIDWN(JZDIM),SIGION(3)
C
      !...... cascade from N2
      IF(JOE.LT.1.OR.JOE.GT.JMAX.OR.E(J).LE.E(JOE)) RETURN
         JOEM1=JOE
         JOEM2=JOE
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-1) JOEM1=JOE+1  !.. May 1999
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-2) JOEM2=JOE+2  !.. May 1999
      CORSIG=SIGEX(3)*EXIFAC
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)
      !.. allocate production from bin J in 2 adjacent bins if small bins
      PREX=0.33333*PREX*CORSIG*XN(3,IZ)
      PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+PREX/DELTE(JOE)
      PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+PREX/DELTE(JOEM1)
      PROCAS(JOEM2,IZ)=PROCAS(JOEM2,IZ)+PREX/DELTE(JOEM2)
      IF(J.EQ.-979) WRITE(6,*) IWR  !.. Dummy statement to use IWR
        RETURN
        END
C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE ACASO2T(JZDIM,IEDIM,J,IZ,IWR,ALT,SIGEX,XN,PROCAS,PHIUP
     > ,PHIDWN,E,DELTE,SIGO1D,SIGION,JMAX,JOE,EXIFAC,BSCX)
C........ to calculate cascade production rates from O2
C........ JOE= index of bin for depositing degraded electrons
      DIMENSION SIGEX(3),XN(3,JZDIM),PROCAS(IEDIM,JZDIM),PHIUP(JZDIM)
     > ,E(IEDIM),DELTE(IEDIM),PHIDWN(JZDIM),SIGION(3)
C
      !...... cascade from O2
      IF(JOE.LT.1.OR.JOE.GT.JMAX.OR.E(J).LE.E(JOE)) RETURN
         JOEM1=JOE
         JOEM2=JOE
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-1) JOEM1=JOE+1  !.. May 1999
         IF(DELTE(JOE).LE.1.0.AND.JOE.LE.JMAX-2) JOEM2=JOE+2  !.. May 1999
      CORSIG=SIGEX(2)*EXIFAC
      PREX=(PHIUP(IZ)*(1-BSCX)+PHIDWN(IZ)*BSCX)*DELTE(J)
      !.. allocate production from bin J in 2 adjacent bins if small bins
      PREX=0.33333*PREX*CORSIG*XN(2,IZ)
      PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+PREX/DELTE(JOE)
      PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+PREX/DELTE(JOEM1)
      PROCAS(JOEM2,IZ)=PROCAS(JOEM2,IZ)+PREX/DELTE(JOEM2)
      IF(J.EQ.-979) WRITE(6,*) IWR  !.. Dummy statement to use IWR
        RETURN
        END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE ACASION(JZDIM,IEDIM,J,IZ,IWR,ALT,XN,PROCAS,PHIUP,PHIDWN
     > ,E,DELTE,SIGION,JMAX,JOE,IONFAC,BSCI)
      !.... calculate secondary and cascade production from ionizing collisions
      REAL IONFAC
      DIMENSION PROCAS(IEDIM,JZDIM),XN(3,JZDIM),SIGION(3),E(IEDIM)
     > ,DELTE(IEDIM)
      !..... Total production of ions at E(J)
      PRION=SIGION(1)*XN(1,IZ)+SIGION(2)*XN(2,IZ)+SIGION(3)*XN(3,IZ)
      PEFLUP=(PHIUP*(1.0-BSCI)+PHIDWN*BSCI)*DELTE(J)*IONFAC
      !..... degraded primaries
      IF(JOE.LE.1.OR.JOE.GE.JMAX) RETURN
      JOEM1=JOE-1
      JOEP1=JOE+1
      IF(JOEM1.GT.J) THEN
         PRION=0.166666666*PRION*PEFLUP
         PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+2.0*PRION/DELTE(JOE)
         PROCAS(JOEM1,IZ)=PROCAS(JOEM1,IZ)+3.0*PRION/DELTE(JOEM1)
         PROCAS(JOEP1,IZ)=PROCAS(JOEP1,IZ)+PRION/DELTE(JOEP1)
      ELSE
         PRION=0.3333333*PRION*PEFLUP
         PROCAS(JOE,IZ)=PROCAS(JOE,IZ)+2.0*PRION/DELTE(JOE)
         PROCAS(JOEP1,IZ)=PROCAS(JOEP1,IZ)+PRION/DELTE(JOEP1)
      ENDIF
      IF(J.EQ.-979) WRITE(6,*) IWR  !.. Dummy statement to use IWR
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C.. calculate secondary production from ionizing collisions. PROBSE=
C.. probability per eV of getting a secondary electron at that energy
C.. Note that the backscatter coefficient (BSCI) determines the
C.. probability of secondary electrons going the opposite direction
C.. to the primary electron. Normally we would assume that the 
C.. angular distribution is isotropic (Opal et al. J. Chem. Phys., 
C.. 1971, page 4100) and therefore BSCI=0.5 for low energy secondaries
      SUBROUTINE ACASSEC(JZDIM,IEDIM,J,IZ,IWR,ALT,XN,PROCAS,PHIUP,PHIDWN
     > ,E,DELTE,SIGION,PROBSE,ISEC,JMAX,BSCI)
      DIMENSION PROCAS(IEDIM,JZDIM),XN(3,JZDIM),SIGION(3),E(IEDIM)
     > ,DELTE(IEDIM),PROBSE(IEDIM)
      !..... PRION= total production of secondary (and primary) ions at E(J)
      PRION=SIGION(1)*XN(1,IZ)+SIGION(2)*XN(2,IZ)+SIGION(3)*XN(3,IZ)
      PEFLUP=(PHIUP*(1.0-BSCI)+PHIDWN*BSCI)*DELTE(J)
      PRION=PRION*PEFLUP
      !... distribute secondaries in appropriate bin
      if (isec.eq.0) isec=isec+1      !xiang changed
      DO 25 IK=ISEC,JMAX
      PROCAS(IK,IZ)=PROCAS(IK,IZ)+PRION*PROBSE(IK)
 25   CONTINUE
      IF(J.EQ.-979) WRITE(6,*) IWR  !.. Dummy statement to use IWR
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C.........calculate the average energy of the secondary electrons from
C........ electrons of energy ELECEN. ELIM= max energy of secondaries,
C........ FNORM= normalization factor to assure total probability =1.
C........  AVESEC= average energy
C******** This routine is no longer used but kept for reference
      SUBROUTINE ASECAVE(ELECEN,AVESEC,EPOTI,SHAPE)
      AVESEC=0.0
      IF(ELECEN.LE.EPOTI) RETURN
         ELIM=0.5*(ELECEN-EPOTI)
         FNORM=1.0/ATAN(ELIM/SHAPE)/SHAPE
         AVESEC=FNORM*0.5*SHAPE**2*ALOG(1.0+(ELIM/SHAPE)**2)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      SUBROUTINE ADIAGS(IEDIM,JMAX,ECHAR,E,EPOTIN,ELPRIM,ELOSSX)
      DIMENSION E(IEDIM),EPOTIN(3,IEDIM),ELOSSX(3,IEDIM)
     > ,ELPRIM(IEDIM)
      !.... average energy loss per ionizing collision
      EREDUC=ECHAR
      DO 126 NCOLL=1,IEDIM
         CALL ASECAVE(EREDUC,AVESEC,20.0,14.0)
         EREDUC=EREDUC-AVESEC-20.0
         IF(EREDUC.LE.30) GO TO 127
C....         WRITE(6,*) NCOLL,EREDUC,AVESEC+20.0
 126  CONTINUE
 127  CONTINUE
      WRITE(117,129)  NCOLL,ECHAR,ECHAR/REAL(NCOLL)
 129  FORMAT(/2X,' # of collisions to thermalize =',I6
     > ,/2X,' Characteristic energy = ',F10.0
     > ,/2X,'Average energy lost per ionizing collision',F9.1)
C
      !.... Printing energy losses per collision
      DO 611 J=1,JMAX
      JC=JMAX+1-J
      WRITE(117,612) E(JC),ELPRIM(JC)-EPOTIN(3,JC)
     >  ,EPOTIN(3,JC),ELPRIM(JC),(ELOSSX(I,JC),I=1,3)
 611  CONTINUE
 612  FORMAT(2X,9F10.2)
      RETURN
      END
C::::::::::::::::::::::::::: PRALTS :::::::::::::::::::::::
       SUBROUTINE ALFLUX(IMAX,IDALT,ALT,INALT,IALT)
       PARAMETER (KDIM=13)
       REAL FLXALT(KDIM),ALT(IDALT)
       INTEGER IALT(KDIM)
       DATA FLXALT/100.,105.,110.,120.,130.,140.,150.,175.,200.
     > ,225.,250.,330.,400./ 
      !...
      K=1
      DO 20 I=2,IMAX-1
         IF(ALT(I-1).LT.FLXALT(K).AND.ALT(I+1).GT.FLXALT(K)) THEN
           IALT(K)=I
           IF (K.LT.KDIM) K=K+1    !xiang changed
         ENDIF
  20   CONTINUE
       INALT=K-1
       !..WRITE(6,*) INALT
       !..WRITE(135,90) (ALT(IALT(I)),I=1,INALT)
C.. 90    FORMAT(5X,' ------- Altitudes for upward fluxes (km) ----------'
C..     > ,/X,'Energy', 22F9.2)
       !..WRITE(136,91) (ALT(IALT(I)),I=1,INALT)
 91    FORMAT(5X,' ------- Altitudes for downward fluxes (km) --------'
     > ,/X,'Energy', 22F9.2)
       RETURN
       END
C:::::::::::::::::::::::::::::: AURFRE :::::::::::::::::::::::::::::
C---- This subroutine calculates the auroral ionization and excitation 
C---- production frequencies for FLIP
      SUBROUTINE AURFRE(IMAX,JMAX,JZDIM,IEDIM,Z,E,PFLIP,ENERAL)
      REAL ISAV,ESAV
      DIMENSION E(IEDIM),PFLIP(IEDIM,JZDIM),XNN(3,401),Z(JZDIM)
     >  ,SIGION(3),SIGEX(3),ENERAL(JZDIM)
      COMMON/AURP2S/ISAV(3,12,401),ESAV(3,12,401),AURZ(401),IAMX


      DO 10 IS=1,3
      DO 10 IK=1,12
      DO 10 IZ=1,401
         ISAV(IS,IK,IZ)=0.0
         ESAV(IS,IK,IZ)=0.0
 10   CONTINUE

      !... densities are set to one to calculate frequency. Will be 
      !... multiplied by density in AURSUB
      DO 15 IZ=1,401
         XNN(1,IZ)=1.0         
         XNN(2,IZ)=1.0         
         XNN(3,IZ)=1.0         
 15   CONTINUE

      !.. Transfer IAMAX and auroral altitude grid for FLIP interpolation
      IAMX=IMAX
      DO I=1,IMAX
        AURZ(I)=Z(I)
        ISAV(1,11,I)=ENERAL(I)  !-- Store energy absorbed by atmosphere
      ENDDO

      !.. LOOP THROUGH ENERGIES CALCULATING PRODUCTION RATES
      DO 20 J=1,JMAX
         IF(E(J).LT.0.8) GO TO 30
         CALL SIGEXS(E(J),SIGEX,SIGION,SIGO1D)
         !..CALL XS_TEST(E(J)) !.. prints cross sections
         !.. LOOP THROUGH FLIP ALTITUDE GRID
         DO 20 I=1,IMAX
           !.. ENERGY LOSS TO COULOMB COLLISIONS WITH THERMAL ELECTRONS
           ET=8.618E-5* 2000.000
           TSIGNE=((3.37E-12)/(E(J)**0.94))*((E(J)-ET)/
     >       (E(J)-(0.53*ET)))**2.36/(E(J)-E(J+1))
           ESAV(1,9,I)=ESAV(1,9,I)+PFLIP(J,I)*TSIGNE
           ESAV(2,12,I)=ESAV(2,12,I)+SIGEX(2)*PFLIP(J,I)
           !.... CALCULATE O PRODUCTION RATES
           CALL ACALOXE(I,E(J),PFLIP(J,I),ESAV)
           !..... CALCULATE N2 EXCITATION RATES
           CALL EPN2XS(I,XNN,PFLIP(J,I),AURZ(I),E(J),SIGEX,SIGION,ESAV)
           !..... CALCULATE ION PRODUCTION RATES
           CALL ACALION(I,E(J),PFLIP(J,I),ISAV)
 20     CONTINUE
 30   CONTINUE


      !---  debug writes
        IF(IMAX.GT.0) RETURN
        DO J=1,IMAX
          WRITE(6,145) J, AURZ(J),ESAV(1,9,J)
     >    ,(ESAV(3,K,J),K=1,4),(ISAV(3,K,J),K=1,3),ISAV(1,11,J)
      ENDDO
 145    FORMAT(1X,I5,F8.1,1P,10E9.1)

      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C...... THIS ROUTINE CALCULATES THE PRODUCTION RATES OF THE O
C...... SPECIES BY SUMMING OVER ENERGY. CROSS SECTIONS ARE FROM OXSIGS
      SUBROUTINE  ACALOXE(I,E,PFLIP,PX)
      DIMENSION SIGEX(22),PX(3,12,401)
      CALL OXSIGS(E,SIGEX,SIGEXT)
      !.. EXCITATION OF O(1D),O(1S), TOTAL-O(1D)
      DO 10 IS=1,6
      PX(1,IS,I)=PX(1,IS,I)+SIGEX(IS)*(PFLIP)
 10   CONTINUE
      PX(1,12,I)=PX(1,12,I)+(SIGEXT-SIGEX(1))*(PFLIP)
      !.. temporary O2 - no partial cross sections
      PX(2,1,I)=PX(2,1,I)+SIGEX(2)*(PFLIP)
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
C...... THIS ROUTINE CALCULATES THE PRODUCTION RATES OF THE IONS
C...... BY SUMMING OVER ENERGY. CROSS SECTIONS FROM SIGEXS
C...... WRITTEN BY P. RICHARDS AUGUST-SEPTEMBER 1988
      SUBROUTINE ACALION(I,E,PFLIP,SECPRD)
      DIMENSION SPRD(3,6),SECPRD(3,12,401),SIGION(3),SIGEXT(3)
      !... BRANCHING RATIOS FOR IONIZATION TO THE VARIOUS STATES OF O+, O2+
      !... AND N2+ BY AURORAL ELECRONS. /4S,O2X,N2X, 2D,A,A, 2P,?,B...
      !... 0.16 IS THE N+ RATIO FROM RICHARDS AND TORR 1985
      !... **** NOTE BRANCHING RATIOS FOR THE HIGHER O+ LEVELS ARE NEEDED
      !--- Branching ratios for ion states were updated Sep 91 by P. Richards
      !--- N2+ from Doering and Goembel, JGR 1991, page 16025. fraction of
      !--- N+ from Richards and Torr JGR 1985, page 9917. fraction of O+ from
      !--- O2 dissociation is from Rapp and Englander Golden J. Chem Phys, 1965
      !--- We still need to find cross sections for the higher levels of O+
      DATA SPRD/.4,.47,.42, .4,.28,.34, .2,.15,.08,0.0,.05,0.0,0.0,
     >   .05,0.0,0.0,0.0,0.16/
C
      !.... IONIZATION RATES FOR O O2, N2. CROSS SECTIONS FROM SIGEXS
      CALL SIGEXS(E,SIGEXT,SIGION,SIGO1D)
 22   DO 130 IS=1,3
      PRION=PFLIP*SIGION(IS)
      DO 25 IK=1,6
      SECPRD(IS,IK,I)=SECPRD(IS,IK,I)+PRION*SPRD(IS,IK)
 25   CONTINUE
 130  CONTINUE
      RETURN
      END
C::::::::::::::::::::::::::::: ATRIS1 ::::::::::::::::::::::::::::::
C
C   Modified by Yingtian Long from Dr. Richards' subroutine
C                    /08/27/97
C
C  This subroutine is used to solve photoelectron flux, phi down, by
C  solving the second order PDE. In th calculation for derivatives
C  at a grid point, the variable spaces between adjacent grid points
C  are considered. The original version of this subroutine is ATRIS0 
C  which can be found in this file.
C
C  Input data PRED is set to zero and BM 1.0 in the calling 
C  subroutine PE2S in the AURORA MODEL program.
C
C... The derivation of the equations was done as a CS 695 course
C... in summer 1997. PGR has a copy of the final paper.

      SUBROUTINE ATRIS1(IDIR,M2,M,J,MT,BM,BG,Z,IMAX,DH,PRED,PRODUP,
     >  PRODWN, PHIDWN,PHIUP,T1,T2,DS)
      PARAMETER (IEDIM=4000,JZDIM = 401)
      REAL DH,Z(JZDIM),BM(JZDIM),BG(JZDIM)
      DIMENSION A(JZDIM),B(JZDIM),C(JZDIM),D(JZDIM),
     >  PRED(JZDIM),PRODWN(IEDIM,JZDIM),PRODUP(IEDIM,JZDIM),
     >  PHIDWN(JZDIM),PHIUP(JZDIM),T1(JZDIM),T2(JZDIM),DS(JZDIM)

      DIMENSION ALPH(JZDIM), GAMM(JZDIM)

      DATA   ISW , FAC2 , JSN , JSW  , KITER , M1 ,ZPRIN ,ZZZZ
     >       /  1  ,   0  ,  3  ,   1   ,  75  ,   1, 300 ,1000 /
      IMXM1=IMAX-1
      !... DO LOOP FOR ITERATING THE SOLUTIONS... PHIDWN IS SOLVED USING
      !... ATRIDAG SOLVER FOR ONE HEMISPHERE, PHIUP IS SOLVED ANALYTICALLY TO UPPER
      !... BDY IN CONJUGATE H-S, PHIUP FOR C.H.S IS FOUND USING ATRIDAG, THEN
      !... PHIDWN IS SOLVED ANALYTICALLY BACK ALONG THE FIELD LINE
      DO 40 I=M2,M
        DELZ=DS(I)+DS(I+1)                            ! ds(i)+ds(i+1)
        DLB=(BM(I+1)-BM(I-1))/(BM(I)*DELZ)            ! (dB)/(B*ds)
        DSLB=2.*((BM(I+1)/DS(I+1)+BM(I-1)/DS(I))/DELZ
     &                  -BM(I)/(DS(I+1)*DS(I)))/BM(I) ! (d^2B)/(B*ds^2)
        DLT1=(T1(I+1)-T1(I-1))/(T1(I)*DELZ)           ! (dT1)/(T1*ds)
        DTS2=(T2(I+1)-T2(I-1))/DELZ                   ! (dT2)/(ds)
        DPR1=(PRED(I+1)-PRED(I-1))/(2*DELZ)           ! 0.5*(dq)/(ds)
        DPR2=(PRODWN(J,I+1)-PRODWN(J,I-1))/DELZ       ! (dq-)/(ds)
        PHI = 1.
        ALPHA = -(DLT1 + 2.*DLB)         ! -[(dT1)/(T1*ds)+2(dB)/(B*ds)]
        BETA = IDIR*(T2(I)*DLT1-DTS2-DSLB)-T2(I)**2+T1(I)**2-ALPHA*DLB
        A(I) = (2.*PHI/DS(I)-ALPHA)/DELZ
        B(I) = BETA-2.*PHI/(DS(I)*DS(I+1))
        C(I) = (2.*PHI/DS(I+1)+ALPHA)/DELZ
        D(I) = (.5*PRED(I)+PRODWN(J,I))*(IDIR*(DLT1+DLB)-T2(I))
     &         -IDIR*(DPR1+DPR2)-T1(I)*(.5*PRED(I)+PRODUP(J,I))
        !!!  where PRED = q, PRODUP = q+ and PRODWN = q-
 40   CONTINUE
C.... END OF D.E. COEFFS --- SOLUTION FOR NEAR H-S ....
      D(M2)=D(M2)-A(M2)*PHIDWN(M2-1)
      D(M)=D(M)-C(M)*PHIDWN(M+1)
      CALL ATRIDAG(JZDIM,PHIDWN,M2,M,A,B,C,D,ALPH, GAMM)

C:::::::: PHIUP IS EVALUATED  ANALYTICALLY   ::::
      DO 60 I=M2,IMXM1
        R1=(T1(I)*PHIDWN(I)+(PRED(I)+2.*PRODUP(J,I))/2.)/T2(I)/BM(I)
        T2DS=T2(I)*DS(I)
        IF(T2DS.GT.70.0)T2DS=70.0
        PHIUP(I)=BM(I)* (R1+(PHIUP(I-1)/BM(I-1)-R1)*EXP(-T2DS))
C....       IF(J.EQ.31.AND.Z(I).GT.250.AND.Z(I).LT.900)
C....     > WRITE(155,93) Z(I),T1(I),T2(I),PRED(I),PRODUP(J,I),BM(I),DS(I)
C....     >  ,R1,PHIDWN(I),PHIUP(I)
C...        IF(I.EQ.IPAS+1) PHIUP(I)=PHIUP(I)*(1-FPAS) ! only needed in FLIP
 60   CONTINUE

      RETURN
 93   FORMAT(1X,F7.0,1P,22E10.1)
      END
C:::::::::::::::::::::::::: APITCH :::::::::::::::::::::::::::::::::::::
C..... Take pitch angle into account assuming adiabatic variation
C..... along the field line. The pitch angle only varies above the upper
C..... boundary altitude for the PDE calculation Z(M). The program could
C..... be written to avoid this if necessary.
      SUBROUTINE APITCH(J,IS,M,M1,IMAX,T1,T2,PRED,PRODUP,PRODWN,BM,Z
     >                 ,DIPANG,AVMU)
      PARAMETER (IEDIM=4000,JZDIM = 401)
      DIMENSION T1(JZDIM),T2(JZDIM),PRED(JZDIM),PRODUP(IEDIM,JZDIM)
     > ,PRODWN(IEDIM,JZDIM), BM(JZDIM),Z(JZDIM)
      !.. RBM used in pitch angle variation
      RBM=0.667/BM(M)
      DO 56 I=M1,IMAX
      !.... pitch angle variation above Z(M)
        IF(Z(I).GT.Z(M)) RAVMU=1.0/SQRT(1.0-BM(I)*RBM)
        IF(Z(I).LT.Z(M)) RAVMU=1.00/SIN(DIPANG)/AVMU
      !... To reset after solution
        IF(IS.EQ.1) RAVMU=1.0/RAVMU
        T1(I)=T1(I)*RAVMU
        T2(I)=T2(I)*RAVMU
        PRED(I)=PRED(I)*RAVMU
        PRODUP(J,I)=PRODUP(J,I)*RAVMU
        PRODWN(J,I)=PRODWN(J,I)*RAVMU
 56   CONTINUE

      RETURN
      END
C::::::::::::::::::::::::::::::::::: EVANS_FLUX ::::::::::::::::::
C...  Auroral fluxes from Evans and Moore, JGR, 1978
      SUBROUTINE EVANS_FLUX(EV,  !.. Evans Energy
     >                    PHIN)  !.. Evans Flux
      INTEGER IADIM
      PARAMETER (IADIM=34)      
      REAL PFLXEV(IADIM), PFLXDWN(IADIM)  !.. Evans energy and down flux
      REAL AVMU         !.. average cosine of pitch angle
      DATA AVMU/0.577/
      DATA PFLXEV/0.5,49.78,64.44,111.08,90.89,135.76,147.95,161.24,
     > 86.08,214.76,262.47,320.79,370.22,465.64,553.02,638.23,780.04,
     > 1038.96,1465.46,2008.65,2127.16,2599.79,3364.89,4112.52,5800.77,
     > 6889.29,7950.88,9176.05,10290.78,10897.94,13319.30,15371.70,
     > 17239.10,18787.06/
      DATA PFLXDWN/5E7,2.135e6,1.118e6,5.853e5,4.861e5,2.751e5,2.123e5,
     > 2.123e5,1.727e5,1.400e5,1.444e5,9.043e4,1.174e5,1.344e5,1.149e5,
     > 1.350e5,1.124e5,1.193e5,1.105e5,1.080e5,8.792e4,7.536e4,6.280e4,
     > 5.112e4,3.561e4,2.606e4,2.010e4,1.363e4,9.734e3,6.782e3,3.812e3,
     > 2.148e3,1.181e3,6.845e2/
      !... get interpolated value of precipitating
      PHIN=0.0
      IF(EV.GE.PFLXEV(1).AND.EV.LE.PFLXEV(IADIM))
     >   PHIN= FTAE(1,IADIM,PFLXEV,PFLXDWN,EV) * AVMU
      RETURN
      END

C::::::::::::::::::::::::::::::::::: FAST_FLUX ::::::::::::::::::
C...  Auroral fluxes from FAST
      SUBROUTINE FAST_FLUX(EV,  !.. FAST Energy
     >                   PHIN)  !.. FAST Flux
      INTEGER IADIM
      PARAMETER (IADIM=48)      
      REAL PFLXEV(IADIM), PFLXDWN(IADIM)  !.. Evans energy and down flux
      REAL AVMU         !.. average cosine of pitch angle
      DATA AVMU/0.577/
      DATA PFLXEV/5.9,9.8,13.7,17.6,21.6,25.5,29.4,35.3,
     > 43.1,51.0,58.8,70.6,86.2,101.9,117.6,141.1,172.5,203.8,
     > 235.2,282.2,345.0,407.7,470.4,564.5,689.9,815.4,940.8,
     > 1129.0,1379.8,1630.7,1881.6,2257.9,2759.7,3261.4,
     > 3763.2,4515.8,5519.4,6522.9,7526.4,9031.7,11038.7,
     > 13045.8,15052.8,18063.4,22077.4,26091.5,30105.6,34119.7/
      DATA PFLXDWN/1.14e8,5.79e7,3.61e7,3.12e7,1.93e7,1.53e7,
     > 1.19e7,7.16e6,4.49e6,3.14e6,2.04e6,1.28e6,8.17e5,7.29e5,6.28e5,
     > 4.82e5,3.30e5,2.98e5,3.10e5,2.54e5,2.26e5,2.16e5,1.71e5,1.66e5,
     > 1.13e5,1.18e5,1.14e5,1.16e5,1.34e5,1.19e5,1.28e5,1.25e5,1.60e5,
     > 1.92e5,3.42e5,5.35e5,2.83e5,9.86e4,2.83e4,8.11e3,2.03e3,7.35e2,
     > 2.12e2,0.00e0,4.75e1,0.00e0,0.00e0,1.59e3/      !... get interpolated value of precipitating
      PHIN=0.0
      IF(EV.GE.PFLXEV(1).AND.EV.LE.PFLXEV(IADIM))
     >   PHIN= FTAE(1,IADIM,PFLXEV,PFLXDWN,EV) * AVMU /6.00
      RETURN
      END
C::::::::::::::::::::::::::::::::::: TJFRE_FLUX ::::::::::::::::::
C...  Auroral fluxes from Fuller-Rowell and Evans JGR, 1987
      SUBROUTINE TJFRE_FLUX(EV,  !.. Evans Energy
     >                    PHIN)  !.. Evans Flux
      INTEGER IADIM
      PARAMETER (IADIM=45)      
      REAL PFLXEV(IADIM), PFLXDWN(IADIM)  !.. Evans energy and down flux
      REAL AVMU         !.. average cosine of pitch angle
      REAL FLXFAC       !.. normalization factor
      DATA AVMU/0.577/
       DATA PFLXEV/385,4.33E+02,4.70E+02,5.25E+02,5.87E+02,6.23E+02,
     > 6.88E+02,7.74E+02,9.46E+02,1.05E+03,1.21E+03,1.39E+03,1.50E+03,
     > 1.68E+03,2.03E+03,2.30E+03,2.53E+03,2.91E+03,3.31E+03,3.70E+03,
     > 4.17E+03,4.91E+03,5.49E+03,6.00E+03,6.71E+03,7.19E+03,7.85E+03,
     > 8.18E+03,8.62E+03,9.09E+03,9.69E+03,1.02E+04,1.09E+04,1.17E+04,
     > 1.21E+04,1.31E+04,1.40E+04,1.52E+04,1.61E+04,1.73E+04,1.91E+04,
     > 2.02E+04,2.15E+04,2.34E+04,2.42E+04/
       DATA PFLXDWN/3.81E+07,3.05E+07,2.56E+07,2.22E+07,1.86E+07,
     > 1.80E+07,1.56E+07,1.38E+07,1.31E+07,1.68E+07,2.37E+07,3.05E+07,
     > 3.68E+07,4.44E+07,5.03E+07,5.88E+07,6.35E+07,6.54E+07,6.13E+07,
     > 5.48E+07,4.75E+07,3.86E+07,3.45E+07,3.13E+07,2.11E+07,1.54E+07,
     > 1.10E+07,8.70E+06,6.55E+06,5.00E+06,3.48E+06,2.66E+06,1.85E+06,
     > 1.30E+06,9.97E+05,6.21E+05,3.93E+05,2.40E+05,1.67E+05,1.07E+05,
     > 6.37E+04,3.91E+04,2.36E+04,1.44E+04,9.57E+03/
      !... get interpolated value of precipitating
      PHIN=0.0
      FLXFAC=.001
      IF(EV.GE.PFLXEV(1).AND.EV.LE.PFLXEV(IADIM))
     >   PHIN= FTAE(1,IADIM,PFLXEV,PFLXDWN,EV) * AVMU * FLXFAC
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::::::::::::: CONSERVE :::::::::::::::::::::::::::
c......This routine calculates and outputs the energy balance for auroral precipitation
C...... P. Richards, June 2019
      SUBROUTINE CONSERVE(M1,    !.. Lower boundary index
     >                  IMAX,    !.. Upper boundary index
     >                 JZDIM,    !.. Energy array dimension
     >                   ZWR,    !.. Altitude for printing
     >                     Z,    !.. Altitude array
     >                 ECHAR,    !.. Characteristic Energy
     >                TENFLX,    !.. Total energy input
     >                  EBSC,    !.. Energy backscattered
     >                  ESUM,    !.. Total precipitating energy
     >                 EHEAT,    !.. Total thermal electron heating
     >                ENABEX,    !.. Total excitation energy  array
     >                ENABIO,    !.. Total ionization energy array
     >                SECPRD)    !.. Ion production rate array
      IMPLICIT NONE
      INTEGER I,IS,M1,IMAX,IEDIM,JZDIM  
      REAL ZWR,Z(JZDIM),ECHAR,EBSC,EHEAT(JZDIM),ENABEX(3,JZDIM),
     >  ENABIO(3,JZDIM),SECPRD(3,JZDIM)
      REAL SUMEH,ELEFT,DELALT,TOTION,EABS,ESUM,ENEBAL,EC
      REAL SUMEX(3),SUMIO(3),TENFLX
      LOGICAL DEBUG
      !... Sum up the energy in the electrons SUMEH, left in the last bin
      !... ELEFT, deposited in excited species (SUMEX), and in ions (SUMIO)
      SUMEH=0.0
      TOTION=0.0
      DO IS=1,3
         SUMEX(IS)=0.0
         SUMIO(IS)=0.0
      ENDDO

      ELEFT=0.0
      DO 120 I=M1+1,IMAX-1
        DELALT=1.0E5*(Z(I+1)-Z(I)) !...+ 0.000*DELZ/BG(I)
        SUMEH=SUMEH+EHEAT(I)*DELALT
        !... Energy left in PRODUP and PRODWN
        !..ELEFT=ELEFT+(PRODUP(JMAX,I)+PRODWN(JMAX,I))*E(JMAX)*DELALT
        DO 140 IS=1,3
          SUMEX(IS)=SUMEX(IS)+ENABEX(IS,I)*DELALT
          SUMIO(IS)=SUMIO(IS)+ENABIO(IS,I)*DELALT
          TOTION=TOTION+SECPRD(IS,I)*DELALT
 140    CONTINUE

        IF(DEBUG) THEN
          IF(I.EQ.M1.OR.MOD(I,15).EQ.0) WRITE(112,202)
 202      FORMAT(2X,'  I     Z     SUMEH    EHEAT    ENABEX   ENABEX',
     >      '   ENABEX   ENABIO   ENABIO   ENABIO')
          IF(DEBUG) WRITE(112,'(1H ,I5,F7.0,1P22E9.1)') I,Z(I),SUMEH
     >      ,EHEAT(I),(ENABEX(IS,I),IS=1,3),(ENABIO(IS,I),IS=1,3)
      ENDIF
 120  CONTINUE

      !.... add up all energy lost to see if it matches the energy deposited
      EABS=SUMEH   !..+ELEFT*0.00
      DO 125 IS=1,3
      EABS=EABS+SUMEX(IS)+SUMIO(IS)
 125  CONTINUE
      IF(ESUM.EQ.0)THEN
         ENEBAL=0.0
      ELSE
         ENEBAL=(EABS+EBSC*0.0)/ESUM
      ENDIF
C
      EC=6.25E11    !*ENEBAL
      WRITE(6,178) ESUM/6.25E11,EABS/EC,EBSC/EC,SUMEH/EC
     >   ,ELEFT/EC,(SUMEX(IS)/EC,IS=1,3),(SUMIO(IS)/EC,IS=1,3)

 178   FORMAT(/2X,' ****** Energy conservation (erg cm-2 s-1)*******'
     > ,/5X,'Energy of precipitating flux                     ',1P,E9.2
     > ,/5X,'Total energy into electrons, ions and neutrals=  ',1P,E9.2
     > ,/5X,'Energy carried by escaping flux=                 ',1P,E9.2
     > ,/5X,'Energy deposited in thermal electrons=           ',1P,E9.2
     > ,/5X,'Energy left in last energy bin=                  ',1P,E9.2
     > ,//5X,'Energy going into excitations of O, O2, N2=     ',1P3E9.2
     > ,/5X,'Energy going into ionizations of O, O2, N2=      ',1P3E9.2)

      WRITE(6,'(/A,A,A/)') 
     >  ' <<<<< NOTE that the escape flux is reflected back down.',
     >  'So, most of the incident energy ends up in the thermosphere.'

      WRITE(6,965) 
     >  SUMEX(1)/SUMIO(1),SUMEX(2)/SUMIO(2),SUMEX(3)/SUMIO(3)
 965  FORMAT(6X,'Ratio going into excitations of O, O2, N2=',1P3E9.2//)

C        IF(ZWR.GT.80) WRITE(6,966) ENEBAL,PHIBSC/PHINSUM
C 966    FORMAT(///,5X,'Energy balance! i.e. Ratio of energy absorbed to'
C     >   ,' energy input=',1P,E9.2,'   RPHIBSC=',1P,E9.2/)


      WRITE(6,967) TOTION/ENEBAL,ESUM/TOTION
 967  FORMAT(/' total # ions= ',1P,E10.2,'  Energy per ion= ',1P,E10.2/)

      WRITE(6,969) 
     > (SUMEX(1)+SUMEX(2)+SUMEX(3))/(SUMIO(1)+SUMIO(2)+SUMIO(3))
 969    FORMAT(/' Ratio SUMEX/SUMIO= ',1P,E10.2/)

      RETURN
      END
C:::::::::::::::::::::::::::::::::::: REESAUR :::::::::::::::::::::::::::::
C....... This routine is based on the Rees, Planet. Space Sci. 1963 paper
      SUBROUTINE REESAUR(ZALT,
     >                     TN,
     >                     OX,
     >                     O2,
     >                     N2,
     >                  ECHAR,
     >                  TEFLX,
     >                PHINSUM,
     >                   RRHO,
     >                   RDEN,
     >                 SUMION,
     >                REESION,
     >                RIONTOT)
      IMPLICIT NONE
      REAL ZALT,TN,OX,O2,N2,ECHAR,TEFLX,PHINSUM,RRHO,RDEN,SUMION,REESION
     >  ,RIONTOT
      REAL ZoverR(101),LAMISO(101),LAMCOS(101),ZMASS,SHGT,Z,R,LAMBDA
      REAL GRAV
       DATA ZoverR/0.00,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,
     > 0.09,0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,
     > 0.20,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.30,
     > 0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.40,0.41,
     > 0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.50,0.51,0.52,
     > 0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.60,0.61,0.62,0.63,
     > 0.64,0.65,0.66,0.67,0.68,0.69,0.70,0.71,0.72,0.73,0.74,
     > 0.75,0.76,0.77,0.78,0.79,0.80,0.81,0.82,0.83,0.84,0.85,
     > 0.86,0.87,0.88,0.89,0.90,0.91,0.92,0.93,0.94,0.95,0.96,
     > 0.97,0.98,0.99,1.00/
       DATA LAMISO/1.32,1.36,1.40,1.44,1.47,1.51,1.54,1.56,
     > 1.58,1.60,1.62,1.63,1.64,1.65,1.66,1.66,1.65,1.64,1.62,
     > 1.60,1.57,1.54,1.51,1.47,1.44,1.40,1.37,1.35,1.32,1.29,
     > 1.26,1.24,1.21,1.19,1.16,1.14,1.12,1.10,1.07,1.05,1.03,
     > 1.00,0.98,0.97,0.96,0.93,0.91,0.88,0.86,0.84,0.81,0.79,
     > 0.77,0.74,0.72,0.69,0.67,0.65,0.62,0.60,0.58,0.55,0.53,
     > 0.51,0.49,0.47,0.45,0.43,0.41,0.39,0.37,0.35,0.34,0.32,
     > 0.30,0.28,0.27,0.25,0.24,0.22,0.21,0.19,0.18,0.16,0.15,
     > 0.14,0.12,0.11,0.10,0.09,0.07,0.06,0.05,0.05,0.04,0.03,
     > 0.02,0.02,0.01,0.005,0.005/
        DATA LAMCOS/1.49,1.49,1.50,1.50,1.50,1.51,1.51,1.51,1.51,
     > 1.51,1.51,1.50,1.50,1.49,1.48,1.47,1.46,1.45,1.45,1.44,1.43,
     > 1.42,1.41,1.40,1.39,1.37,1.36,1.34,1.33,1.31,1.30,1.28,1.27,
     > 1.25,1.23,1.21,1.19,1.17,1.15,1.13,1.11,1.09,1.07,1.05,1.03,
     > 1.01,0.98,0.96,0.94,0.91,0.90,0.89,0.87,0.84,0.82,0.79,0.77,
     > 0.74,0.72,0.69,0.67,0.64,0.62,0.60,0.57,0.55,0.52,0.50,0.48,
     > 0.45,0.43,0.41,0.39,0.37,0.35,0.33,0.31,0.30,0.28,0.26,0.24,
     > 0.23,0.21,0.19,0.18,0.16,0.15,0.13,0.12,0.11,0.10,0.09,0.08,
     > 0.07,0.06,0.05,0.04,0.04,0.03,0.02,0.00/

      ZMASS=1.67E-24*(16*OX+32*O2+28*N2)/(OX+O2+N2)
      GRAV=980.665/(1.0+ZALT/6370.0)**2  !.. gravity
      SHGT=1.38E-16*TN/GRAV/ZMASS
      Z=(OX+O2+N2)*ZMASS*SHGT
      R=4.57E-6*(ECHAR/1000)**1.75
      LAMBDA=1.32
      IF(Z/R.GT.0.01.AND.Z/R.LT.1.0) THEN
        LAMBDA=LAMISO(NINT(100*Z/R)+1)
      ELSEIF(Z/R.GE.1.0) THEN
        LAMBDA=0.00
      ENDIF
      REESION=(ECHAR/35)*(RRHO/R)*LAMBDA*(OX+O2+N2)/RDEN
      WRITE(6,'(F10.2,1P,22E10.2)') 
     > ZALT,RRHO,ZMASS*(OX+O2+N2),(OX+O2+N2),Z,R,Z/R,LAMBDA,
     > REESION,REESION*PHINSUM,SUMION,PHINSUM,
     > RIONTOT/(SUMION+.00001),RIONTOT

      RETURN
      END
C::::::::::::::::::::::::::: FITREES ::::::::::::::::::::::::::::::::::
C.......exponential interpolation of production rates (fitting y= A exp(Bx)
      REAL FUNCTION FITREES(IDIM,JMIN0,JMAX0,ZP,T,Z)
      IMPLICIT NONE
      INTEGER K,IDIM,JMIN0,JMAX0
      REAL ZP(IDIM),T(IDIM),Z,TK,TK1,B

      !.. Use bisection to find the nearest altitude ZP(K) to desired 
      !.. altitude (Z)
       CALL BISPLIT(JMIN0,JMAX0,IDIM,Z,ZP,K)

      !.. Now interpolate to get desired production rate. Test to make
      !.. sure we haven't run off the end of the field line grid
 2    IF(K.EQ.JMAX0) THEN
         FITREES=T(K)
      ELSE IF(T(K).LT.0.OR.T(K+1).LT.0) THEN
         FITREES=0.0
      ELSE 
         TK=(T(K)) + 1.0E-22
         TK1=(T(K+1)) + 1.0E-22
         B=ALOG(TK/TK1)/(ZP(K)-ZP(K+1))
         FITREES=TK*EXP(B*(Z-ZP(K)))
         IF(K.EQ.1.AND.Z.EQ.ZP(K)) FITREES=TK
      ENDIF
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::: TIROGAUSS :::::::::::::::::::::
C...... This routine simulates Tim's TIROS data using narrow Gaussians.
C...... The results can be compared to ionization rates using the Rees
C...... formulation
      SUBROUTINE TIROGAUSS(EV,TENFLUX,PHIN)
      IMPLICIT NONE
      INTEGER I
      REAL EV,TENFLUX,PHIN,ECHAR,WEQ,SFAC,TENFLX
      REAL ETIRO(12),PHITIRO(12),PHINORM
      DATA ETIRO/24469,17293,12141,8621,6043,4198,2917,2047,1394,938,
     >   604,384/
c      DATA PHITIRO/0.0002,0.0028,0.0240,0.18,0.80,1.13,1.68,1.26,0.78,
c     >   0.33,0.5,1.0/
c      DATA PHINORM/7.68/   !.. normalizes PHITIRO
      DATA PHITIRO/0.02,0.56,2.4,7.2,16,14,10.08,3.78,0.78,0.33,
     > 0.25,0.2/
      DATA PHINORM/55.6/   !.. normalizes PHITIRO
     
      DATA SFAC/0.05/
      PHIN=0.0
      DO I=1,12
        ECHAR=ETIRO(I)
        TENFLX=1.00 * TENFLUX*PHITIRO(I)/PHINORM
        WEQ=SFAC*ECHAR
        !..WEQ=30.0
        PHIN=PHIN+0.564*TENFLX*EXP(-(EV-ECHAR)**2/WEQ**2)/(WEQ*ECHAR)
      ENDDO
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::: TIROHIST :::::::::::::::::::::
C...... This routine simulates Tim's TIROS data using Histograms.
C...... TJFR email 2019-07-19 The energy ranges for the 15 bins are:
C....... Data EN/.37,.6,.92,1.37,2.01,2.91,4.19,6.,8.56,12.18,
C....... 17.3,24.49,36.66,54.77,81.82/
C....... data WIDTH/.158,.315,.315,.63,.631,1.261,1.26,2.522,
C....... 2.522,5.043,5.043,10.,14.81,22.13,33.06/
C....... EN is the mean energy in each bin, and WIDTH is the width of the bin.
C....... So EN x WIDTH is the energy influx, normalized to 1 erg/cn2/s. So is 
C....... scaled by whatever the local total energy at that location from the 
C....... TIROS energy influx maps. To get the boundaries of the bins start 
C....... from 0.3 and add the width successively.
C....... The results can be compared to ionization rates using the Rees
C....... formulation
      SUBROUTINE TIROHIST(ITIRO,EV,AVMU,PHIN)
      IMPLICIT NONE
      INTEGER I,IDIM,ITIRO
      PARAMETER (IDIM=13)
      REAL EV,AVMU,PHIN,ECHAR
      REAL EBDY(IDIM),ETIRO(IDIM),PHITIRO(IDIM),PHINORM(IDIM)
      DATA EBDY/300,458,773,1088,1718,2349,3610,4870,7392,9914,14957,
     >   20000,30000/
      DATA PHITIRO/1.240E+05,6.183E+04,4.090E+04,9.630E+04,1.557E+05,
     > 2.086E+05,1.403E+05,9.893E+04,2.189E+04,2.980E+03,3.470E+02,
     > 30.0,0.0/
      DATA ETIRO/379,616,931,1403,2034,2980,4240,6131,8653,12436,
     >   17479,25000,30000/
      DATA PHINORM/13*0.285/   !.. normalizes PHITIRO

      !.. To select one energy bin, comment the loop and set ITIRO
      !.. ITIRO=1
      !.. I=ITIRO 
      PHIN=0.0
      DO I=1,IDIM-1     
        IF(EV.GE.EBDY(I).AND.EV.LT.EBDY(I+1)) THEN
          !... AVMU here because PHIN is divided by it in main
          PHIN=PHITIRO(I) * PHINORM(I) * AVMU
          ECHAR=ETIRO(I)
        ENDIF
      ENDDO      
      RETURN
      END
C::::::::::::::::::::::::::::::::::::::: PHINREES :::::::::::::::::::::
C...... This routine puts the 2-stream fluxes into the bins for the REES
C...... procedure
C...... TJFR email 2019-07-19 The energy ranges for the 15 bins are:
C....... Data EN/.37,.6,.92,1.37,2.01,2.91,4.19,6.,8.56,12.18,
C....... 17.3,24.49,36.66,54.77,81.82/
C....... data WIDTH/.158,.315,.315,.63,.631,1.261,1.26,2.522,
C....... 2.522,5.043,5.043,10.,14.81,22.13,33.06/
C....... EN is the mean energy in each bin, and WIDTH is the width of the bin.
C....... So EN x WIDTH is the energy influx, normalized to 1 erg/cn2/s. So is 
C....... scaled by whatever the local total energy at that location from the 
C....... TIROS energy influx maps. To get the boundaries of the bins start 
C....... from 0.3 and add the width successively.
C....... The results can be compared to ionization rates using the Rees
C....... formulation
      SUBROUTINE PHINREES(EV,DELTE,PHIN,PHITIRO,ETIROS)
      IMPLICIT NONE
      INTEGER I,IDIM
      PARAMETER(IDIM=13)
      REAL EV,DELTE,PHIN,PHITIRO(IDIM),ETIROS(IDIM)
      REAL EBDY(IDIM),ETIRO(IDIM)
      DATA EBDY/300,458,773,1088,1718,2349,3610,4870,7392,9914,14957,
     >   20000,30000/
      DATA ETIRO/379,616,931,1403,2034,2980,4240,6131,8653,12436,
     >   17479,25000,30000/
     
      !... Add FLIP flux to correct energy interval
      DO I=1,IDIM-1     
        IF(EV.GE.EBDY(I).AND.EV.LT.EBDY(I+1)) THEN
          PHITIRO(I)=PHITIRO(I)+PHIN*DELTE
          ETIROS(I)=ETIRO(I)
        ENDIF
      ENDDO      
      RETURN
      END
      
 
      
 
